home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / toolbox / czt.r < prev    next >
Text File  |  1994-04-25  |  2KB  |  92 lines
  1. //-------------------------------------------------------------------//
  2. // czt.r
  3.  
  4. // Syntax:    czt ( X , M , W , A )
  5.  
  6. // Description:
  7.  
  8. //    czt returns the M element Chirp z-transform, where:
  9. //    M: The number of points in the returned spectrum.
  10. //    W: The "spacing" between points along the z-pane spiral.
  11. //       A complex number along the unit circle.
  12. //    A: the complex starting point of the evaluated transform.
  13. //       A complex number along the unit circle.
  14.  
  15. //    At present X must be a row or column matrix. The column-wise
  16. //    matrix czt operation will be added later.
  17.  
  18. //      Ian Searle, 10/10/93
  19.  
  20. //    References:
  21. //    Rabiner, L.R. and B. Gold, Theory and Application of
  22. //     Digital Signal Processing, Prentice-Hall, Englewood
  23. //     Cliffs, New Jersey, pp. 393-399, 1975.
  24. //-------------------------------------------------------------------//
  25.  
  26. czt = function ( _x, M, W, A )
  27. {
  28.   local (G, K, L, N, V, X, Y, ...
  29.          column, g, k, n, nr, nc, v, x, y, y1, y2)
  30.  
  31.   x = _x;
  32.   if (min (size (x)) != 1) { error ("czt: only vector X allowed"); }
  33.   nr = x.nr; nc = x.nc;
  34.  
  35.   if (nc == 1) 
  36.   { 
  37.     x = x'; 
  38.     nr = x.nr; nc = x.nc;
  39.     column = 1;
  40.   else
  41.     column = 0;
  42.   }
  43.   N = nc;
  44.  
  45.   if (!exist (M)) { M = length (x); }
  46.   if (!exist (W)) { W = exp (-1j .* 2 .* pi ./ M); }
  47.   if (!exist (A)) { A = 1; }
  48.  
  49.   //
  50.   // Compute a power of 2 value for fft()
  51.   //
  52.  
  53.   L = 1;
  54.   while (L < (N + M - 1)) { L = 2 .* L; }
  55.  
  56.   //
  57.   // Form the L-point sequence y(n)
  58.   //
  59.  
  60.   n = 0:(N-1);
  61.   y = [ A.^(-n) .* W.^((n.^2)./2) .* x , zeros (size (N+1:L)) ];
  62.  
  63.   //
  64.   // Form the L-point sequence v(n)
  65.   //
  66.  
  67.   n = 0:(M-1);
  68.   v = W.^(-(n.^2)./2);
  69.   v = [v, zeros (size (M+1:(L - N + 1)))];
  70.   n = (L-N+1):(L-1);
  71.   v[(L-N+2):L] = W.^(-((L-n).^2)./2);
  72.  
  73.   //
  74.   // Convolution
  75.   //
  76.  
  77.   Y = fft (y, L);
  78.   V = fft (v, L);
  79.   G = V.*Y;
  80.   g = ifft (G, L);
  81.  
  82.   k = 0:(M-1);
  83.   X = g[1:M] .* W.^((k.^2)./2);
  84.  
  85.   if (column)
  86.   {
  87.     return X';
  88.   else
  89.     return X;
  90.   }
  91. };
  92.